function MyTimeFrequency
%MYTIMEFREQUENCY Generate timefrequency data files
%  Uses information in Database.timefrequencyrequests and
%  writes result to Database.timefrequency.



%% Initialize
global Database;

%% If timefrequencyrequest doesn't exist, create it..
if ~any(strcmp(fieldnames(Database),'timefrequencyrequests'))
    frequencyRange = 4:200 %#ok<NOPRT>
    timesOut = 200 %#ok<NOPRT>
    GenerateTimeFrequencyRequestsTable(frequencyRange, timesOut);
end

%% Timefreq the requested items

cSubject = strcmp(Database.timefrequencyrequests_header, 'Subject');
cCondition = strcmp(Database.timefrequencyrequests_header, 'Condition');
cChannel = strcmp(Database.timefrequencyrequests_header, 'Channel');
cChannelType = strcmp(Database.timefrequencyrequests_header, 'ChannelType');
cSourcefile = strcmp(Database.timefrequencyrequests_header, 'Sourcefile');
cFrequencyRange = strcmp(Database.timefrequencyrequests_header, 'FrequencyRange');
cTimesOut = strcmp(Database.timefrequencyrequests_header, 'TimesOut');
% Select the distinct sourcefiles
sourcefiles = unique(Database.timefrequencyrequests(:,cSourcefile)); % TODO: Unique makes that perhaps not all request will be processed..
for iS = 1:size(sourcefiles)
    % Load the file
    EEG = pop_loadset('filename', sourcefiles{iS});
    % Get all the channels that need to be processed for this one dataset..
    rSourcefile = strcmp(Database.timefrequencyrequests(:,cSourcefile), sourcefiles{iS});
    for i = find(rSourcefile)'
        % Get all the variables
        subject = Database.timefrequencyrequests{i,cSubject};
        condition = Database.timefrequencyrequests{i,cCondition};
        channel = Database.timefrequencyrequests{i,cChannel};
        channelType = Database.timefrequencyrequests{i,cChannelType};
        frequencyrange = Database.timefrequencyrequests{i,cFrequencyRange};
        timesout = Database.timefrequencyrequests{i,cTimesOut};
        % Create the filename
        filename = [Database.Path.TimeFrequencyData subject ' ' channel ' ' condition ' ' num2str(min(frequencyrange)) '-' num2str(max(frequencyrange)) 'x' num2str(timesout) '.mat']; %#ok<NOPRT>
        disp([filename '...']);
        % Select the data
        selectedElectrode = strcmp({EEG.chanlocs.labels}, channel);
        if sum(selectedElectrode)==1 % There can only be one!
            data = squeeze(EEG.data(selectedElectrode,:,:));
            samplingrate = EEG.srate;
            if isempty(EEG.times), timesAdjust = EEG.xmin * 1000; else timesAdjust = min(EEG.times); end % Just in case there is only 1 epoch (again) (Note:?????)
            [freqs times nrOfEpochs] = TimeFreqThisToFile( data, samplingrate, timesAdjust,  ...
                subject, condition, channel, channelType, frequencyrange, timesout, ...
                filename);
            UpdateOrAddTimeFrequencyFileInDatabase(filename, subject, channel, channelType, condition, freqs, times, nrOfEpochs)
        else
            warning ('File skipped!'); %#ok<WNTAG>
        end
    end
end



end

%% TF this channel
function [frequencies times nrOfEpochs] = TimeFreqThisToFile(eegdata, samplingrate, timesAdjust, subject, condition, channel, channelType, frequencyRange, timesout, filename)

[data frequencies times coherence] = timefreq(eegdata, 512, ...
    'ntimesout', timesout, ...
    'freqs', frequencyRange, ...
    'cycles', [3 0.5], ...
    'wletmethod', 'dftfilt3');
power = data.*conj(data); % power for wavelets
times = times + timesAdjust; % Correction for t=0, since timefreq fucks up the offset..

%% Normalization
power = NormalizePower(power, times);  
disp([subject ' - ' condition ' - ' channel ': TF power normalized..']);
nrOfEpochs = size(power,3); %#ok<NASGU>

save(filename, 'power', 'coherence', 'frequencies', 'times', 'channel', 'channelType', 'subject', 'condition', 'nrOfEpochs', '-v7.3');
end


% [alltfX freqs timesout R] = timefreq(data, g.srate, tmioutopt{:}, ...
%     'winsize', g.winsize, 'tlimits', g.tlimits, 'detrend', g.detrend, ...
%     'itctype', g.type, 'subitc', g.subitc, 'wavelet', g.cycles, ...
%     'padratio', g.padratio, 'freqs', g.freqs, 'freqscale', g.freqscale, ...
%     'nfreqs', g.nfreqs, 'timestretch', {g.timeStretchMarks', g.timeStretchRefs}, 'wletmethod', g.wletmethod);



%% Generate a new table for tf requests
function GenerateTimeFrequencyRequestsTable(frequencyRange, timesOut)

global Database;

Database.timefrequencyrequests_description = ...
    'Contains requests (channels, ranges, times) to run for tf processing';
Database.timefrequencyrequests_header = ...
    {'Subject', 'Condition', 'Channel', 'ChannelType', 'Sourcefile', 'FrequencyRange', 'TimesOut'};
cSubject = strcmp(Database.timefrequencyrequests_header, 'Subject');
cCondition = strcmp(Database.timefrequencyrequests_header, 'Condition');
cSourcefile = strcmp(Database.timefrequencyrequests_header, 'Sourcefile');
cChannel = strcmp(Database.timefrequencyrequests_header, 'Channel');
cChannelType = strcmp(Database.timefrequencyrequests_header, 'ChannelType');
cFrequencyRange = strcmp(Database.timefrequencyrequests_header, 'FrequencyRange');
cTimesOut = strcmp(Database.timefrequencyrequests_header, 'TimesOut');

includeRegions = GetRegions();
includeSubjects = {'CW', 'DS', 'JF'};
dt = cell(0, size(Database.timefrequencyrequests_header, 2));
for iS = 1:length(includeSubjects)
    subject = includeSubjects{iS};
    unsortedElecs = {};
    % Find all the distinct electrodes
    for iR = 1:length(includeRegions)
        region = includeRegions{iR};
        unsortedElecs = [unsortedElecs GetRegionElectrodes(subject, region)];
    end
    electrodes = unique(unsortedElecs);
    
    % Find all the conditions for the subject
    cEpochsSubject = strcmp(Database.epochs_header, 'Subject');
    cEpochsCondition = strcmp(Database.epochs_header, 'Short');
    cEpochsFilename = strcmp(Database.epochs_header, 'Filename');
    rEpochsSelected = strcmp(Database.epochs(:,cEpochsSubject), subject);
    conditions = Database.epochs(rEpochsSelected, cEpochsCondition);
    sourcefiles = Database.epochs(rEpochsSelected, cEpochsFilename);
    
    % Now add permutation to ToDo-datatable
    for iCo = 1:length(conditions)
        for iEl = 1:length(electrodes)
            newRow = cell(size(Database.timefrequencyrequests_header));
            newRow{1,cSubject} = subject;
            newRow{1,cCondition} = conditions{iCo};
            newRow{1,cSourcefile} = sourcefiles{iCo};
            newRow{1,cChannel} = electrodes{iEl};
            newRow{1,cChannelType} = 'Electrode';
            newRow{1,cFrequencyRange} = frequencyRange;
            newRow{1,cTimesOut} = timesOut;
            dt = [dt;newRow]; %#ok<AGROW>
        end
    end
end
Database.timefrequencyrequests = dt;
end